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Abstract 

We propose a new method for molecular dynamics and Monte Carlo simulations, which is referred 
to as the replica-permutation method (RPM), to realize more efficient sampling than the replica- 
exchange method (REM). In RPM not only exchanges between two replicas but also permutations 
among more than two replicas are performed. Furthermore, instead of the Metropolis algorithm, 
the Suwa-Todo algorithm is employed for replica-permutation trials to minimize its rejection ratio. 
We applied RPM to particles in a double-well potential energy, Met-enkephalin in vacuum, and a 
C-peptide analog of ribonuclease A in explicit water. For a comparison purposes, replica-exchange 
molecular dynamics simulations were also performed. As a result, RPM sampled not only the 
temperature space but also the conformational space more efficiently than REM for all systems. 
From our simulations of C-peptide, we obtained the a-helix structure with salt-bridges between 
Gly2 and ArglO which is known in experiments. Calculating its free-energy landscape, the folding 
pathway was revealed from an extended structure to the a-helix structure with the salt-bridges. 
We found that the folding pathway consists of the two steps: The first step is the "salt-bridge 
formation step" , and the second step is the "a-helix formation step" . 
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I. INTRODUCTION 



In recent years, generalized-ensemble algorithms are frequently employed to study 
biomolecular systems (for reviews, see, e.g., Refs. This is because it is difficult 

to obtain sufficient sampling in the conformational space of these systems by conventional 
canonical-ensemble simulations 3|-|9|]. The canonical simulations tend to get trapped in a 
local-minimum free-energy state of the biomolecular systems. 
The replica-exchange method (REM) 



IG 



11| is one of the most well-known methods 



among the generalized-ensemble algorithms. Non-interacting copies (replicas) of a target 
system are employed in REM. Different temperatures are assigned to the replicas. By ex- 
changing the temperatures between the replicas, random walks of the replicas in the temper- 
ature space are realized. Accordingly, the simulation can escape from local-minimum states. 



It is easier to perform replica-exchange molecu 
ulations than to perform multicanonical MD 



ar dynamics (MD) or Monte Carlo (MC) sim 
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13| or MC simulations 



although 



the multicanonical algorithm is also one of the most wel 



-known generalized-ensemble algo- 



16|-|25| non-Boltzmann weight factors 



rithms. In the multicanonical and similar algorithms 
are used as the weight factors. These non-Boltzmann weight factors are not a priori known 
and have to be determined by iterative procedures. In REM, however, the usual Boltzmann 
weight factor is employed for each replica. Therefore, there is no necessity to perform the 
procedures of obtaining the weight factor. 

In REM the Metropolis algorithm is utilized to exchange the temperatures between 
the replicas. The Metropolis algorithm is a Markov chain Monte Carlo (MCMC) method 
and widely employed to obtain a required ensemble. The Metropolis algorithm is designed 
so as to satisfy the detailed balance condition, which is a sufficient condition to perform 
state transitions based on MCMC. Trials of such state transitions are accepted or rejected 
stochastically, and their rejection ratio increases when there is a large difference between 
the probability distributions. In order to minimize the rejection ratio, new algorithm was 
proposed recently by Suwa and Todo |26 |. 

In the Suwa- Todo algorithm the detailed balance condition is not imposed. They intro- 
duced a graphical procedure called weight allocation instead of solving the detailed balance 
condition algebraically. By minimizing the rejection ratio with the weight allocation, this al- 
gorithm realizes efficient sampling of states. For a system which has only two states for each 
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particle like the Ising model, however, this algorithm is exactly the same as the Metropolis 



algorithm 



26|. 



To take advantage of the Suwa-Todo algorithm, one might consider to exchange 
temperatures in REM by this algorithm. For example, let us consider the follow- 
ing exchange of the temperatures Ti and T2 between the replicas. Replica 1 and 
Replica 2: (Replica 1 at Ti; Replica 2 at T2) — )■ (Replica 1 at T2; Rephca 2 at Ti). In this 
case, the number of all combinations of the replicas and temperatures is only two, 
(Replica 1 at Ti; Replica 2 at T2) and (Replica 1 at T2; Replica 2 at Ti). For this exchange, 
therefore, the Suwa-Todo algorithm is exactly the same as the Metropolis algorithm. In 
general, because temperatures are exchanged between two replicas as in this example, REM 
is not able to take advantage of the Suwa-Todo algorithm. 

In order to solve this difficulty, we propose a new generalize-ensemble algorithm, which 
we refer to as the replica-permutation method (RPM). In this method temperatures are 
permuted among more than two replicas by the Suwa-Todo algorithm. We can utilize the 
Suwa-Todo algorithm because the number of all combinations of the replicas and tempera- 
tures is larger than two. 

We first apply this new method to particles in a double-well potential energy. For a 
comparison purposes, RPM with the Metropolis algorithm, and REM are also performed. As 
the second application of RPM, we employ Met-enkephalin in vacuum. This penta-peptide 
is often used as a test system to see the usefulness of new algorithms 27H30|. Furthermore, 
RPM is applied to a C-peptide analog of ribonuclease A in explicit water, which is known 



to form a a-helix structure 



31 



lJ-|35|], to see its sampling efficiency for a larger biomolecular 
system. The results of the second and third applications are compared with those of REM. 

It is considered that the a-helix structure of the C-peptide analog is stabilized by salt 
bridges (SBs) between Gly2 and ArglO 3J]. We discuss the role of the SBs for the a-helix 
structure. Furthermore, a folding pathway based on a free-energy landscape is presented 
from our simulations. 

In Section [Til we describe formulation of RPM. The Suwa-Todo algorithm and the graph- 
ical procedure, weight allocation, for MCMC are also introduced in this section. We present 
the details of our simulations in Section lllli The results are shown in Section [IVl Section |V] 
is devoted to conclusions. 
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II. METHODS 



A. Markov Chain Monte Carlo Method with the Weight Allocation 

We first describe usual MCMC method with the Metropolis and Suwa-Todo algorithms. 
Let us consider that a system has n states and that the weight of state i is expressed by Wi 
(i = 1, ■ ■ ■ , n). In MCMC the weight Wi is given by 

n 

Wi = ^WjP{] ^i) , (1) 

i=i 

where P{j — )■ i) is the transition probability from state j to state i. By defining the amount 
of stochastic flow from state j to state i as 

v{] ^i)= WjP{j -> i) , (2) 

Eq. ([T]) can be rewritten as 

n 

w, = ^v{j i) . (3) 

i=i 

This amount of stochastic flow v{j — )■ i) also satisfies 

n 
1=1 

because 

n 

5^P(j^z) = l. (5) 

j=i 

From Eqs. (E]) and (jl]), the global balance equation is derived: 

n n 

J^v^i^ j) = J2v^j ^i) , (6) 

1=1 i=l 

By performing state transitions with v{j — > i) satisfying this equation, a required ensemble 
is obtained. 

In the Metropolis algorithm the detailed balance condition, which is a sufficient condition 
for Eq. (|6]), is employed to obtain the required ensemble. Here, the detailed balance condition 
is described in terms of the amount of stochastic flow as follows: 

vii j) = v{j i) . (7) 
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v{j — )• i) is given by the following equation so as to satisfy Eq. ([7]): 

-^i) = ^^-^min [wj, Wi] , j 7^ i , (8) 

where the coefficient l/{n — 1) comes from a random selection of state i from {n — 1) 
candidates except state j. 

When v{j — )• i) satisfies Eqs. (j3]) and (jlj), Eq. is automatically fulfilled . Therefore, 
we focus on Eqs. ([3]) and (jl]) after this. These two equations can be understood visually 
by the weight allocation as follows (see Fig. [U^a)): Equation ([3]) is regarded as filling a box 
whose size is Wi by blocks whose sizes are v{j — )■ i) (j = I,-- - ,n) without any vacant 
space. Equation (j4]) is regarded as dividing a block whose size is Wj into smaller blocks 
whose sizes are v{j i) {i = 1, ■ ■ ■ ,n). In order to satisfy both equations simultaneously 
for any i and j, therefore, we prepare blocks Wj and boxes Wi {i,j = 1, ■ ■ ■ ,n). By using 
blocks v{j — 7- i) divided from wj, all Wi are filled without any space. Figure [l](a) shows this 
weight allocation of Metropolis algorithm for a system which has four states (n = 4). The 
block size of v{j i) is calculated from Eq. (jH]). Red frame blocks represent rejected fiow 
f (i — )■ i) (i = 1, ■ ■ ■ , 4), and their sizes are associated directly with the average rejection 

ratio Y.i'"{^-^ ^)/Y.i^i- 

A new algorithm is proposed recently by Suwa and Todo through the weight allocation 
to minimize the average rejection ratio for state transitions in MCMC We refer to this 
algorithm as the Suwa- Todo algorithm. In this algorithm v{j — t- i) satisfies Eqs. (JSj) and 
(jl]) without imposing the detailed balance condition in Eq. ([7]). In the Suwa- Todo weight 
allocation, wj is divided and Wi is filled as follows (see Fig. [^b)): 

(i) Choose the state which has the maximum weight. If two or more states have the 
maximum weight, one of them is chosen. Here, we assume that Wi is the maximum 
weight without loss of generality. 

(ii) Box W2 is filled by block wi (f (1 — )■ 2)). If block wi still remains after filling box W2, 
try to fill the next box W3 {v{l — 3)). this process is continued until the block size of 
wi become {v{l — > 4), ■ ■ ■ , v{l — )■ k)). 

(iii) By using block W2, fill boxes in turns from the last partially filled box at Step (ii) 
(f(2 — !■ k),--- ,t>(2 — i- /)). This procedure is repeated for the blocks W3, ■ ■ ■ 

(^(3->0,---)- 
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(iv) Once all the boxes except wi are saturated, box wi (■ ■ ■ ,v{n — )■ 1)) is filled. 

Figure [U^b) shows the weight allocation for n = 4. In this figure the average rejection 
ratio (z — ?► 0- ■'■^ ^^^^ Suwa-Todo algorithm, v{j — )■ i) is given by 

v{j —> i) = maLx[0,mm[Aji,Wj + Wi — Aji,Wj,Wi]] , (9) 



where 



and 



Aji = Sj - S',_i + wi , (10) 



S*, = ^ Wj , 5*0 = . (11) 

i=i 

By performing state transitions based on the amount of stochastic flow v{j — ?■ i) or the 
transition probability P{j i) = v{j — i- i)/wj, the required ensemble is obtained. 
Regarding the rejection ratio, from Eq. ([9]), 

I meLx[0,2wi- Sn] , i = I , 
v[i (12) 

[ , i>2 . 

Therefore, the rejection ratio becomes 0, if 

^1 < Y ■ (13) 

It means a reject-free MC simulation is realized if Wi is less than or equal to the half of the 
total weight Sn = X]j=i "^j ■ This condition can be fulfilled in most cases. 

B. Replica-Exchange Method with the Metropolis Algorithm 

We consider a system of atoms with their coordinate vectors and momentum vectors 
denoted by g = {q^, ■ ■ ■ , q^} and p = {pi, ■ ■ ■ ,Pn}, respectively. The Hamiltonian H in 
state X = {q,p) is given by the sum of the kinetic energy K and potential energy V: 

H{x) = K{p) + V{q) . (14) 

In the canonical ensemble at temperature T, each state x is weighted by the Boltzmann 
factor: 

W^{x) = e-^^(") , (15) 



where P = l/fce^ (^b is the Boltzmann constant). 

Let us suppose that there are M non-interacting copies (or rephcas) of the original system 
in the canonical ensemble at M different temperatures {m = 1, ■ ■ ■ ,M). In REM the 
replicas are arranged so that there would be always exactly one replica at each temperature. 
In other words, there is a one-to-one correspondence between the replicas and temperatures. 
Therefore, the label i {i = 1, ■ ■ ■ , M) for the replicas is a permutation of the label m 
(m = 1, ■ ■ ■ , M) for the temperatures, and vice versa: 

i = i(m) = f(m) , 
m = m{i) = f~^{i) , 

where f{m) is a permutation function of m and f^^{i) is its inverse. 

Let = jxp^ ■ ■ ■ = \^x^^(^^y ■ ■ ■ ,x^^[[^^'j stands for a "state" in REM. Here, 

the superscript i and the subscript m in Xm are labels of the replicas and temperatures, 
respectively. All possible combinations between the replicas and temperatures are labeled 
by the subscript a. The state is specified by the M sets of coordinates g'*' and momenta 

of atoms in replica i at temperature Tm- 



XI 



il^^P^'^)^ ■ (17) 



Because the replicas are non-interacting, the weight factor wr(Xq,) for the state Xa is given 
by the product of Boltzmann factors for each replica i (or at temperature Tm)'- 

M M 

(X„) = JJexp {-Pm(i)H } = H^^P (-/^^^ (^1^^"^^') } ' (18) 

i=l m=l 

where i{m) and m{i) are the permutation functions in Eq. (fT6i) . 

We now consider exchanging a pair of replicas in REM. Suppose we exchange replicas j 
and k which are at temperatures Tm and T„, respectively (j = z(m), k = i{n)): 

Xq, { ■ ■ ■ ? -^m ; ■ ■ ■ ; -^L ^ ? ' ' ' } ^ { ' ' ' ' "^m ? ' ' ' ; -^n^ ' ' ' ' } ' (1^) 

From Eq. ([8]) the amount of stochastic flow v {X^ — )■ X^) for this replica exchange is given 
by 

V (X, -> X^) = Cmin [wn (X,) , (X^)] , (20) 
and the transition probability P (Xq, -t- X^) is expressed by 



P (X„ ^ Xp) = Cmm 



(21) 
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where C = 1/ (mC2) if replicas j and k were selected randomly among M replicas. Here, 
A/C2 is the number of 2-combinations from M elements: = M!/{(M — 2)!2!} . In the 
usual REM, replica exchanges are tried only between neighboring two temperatures. In this 
case, C = 1 because M = 2. 



C. Replica- Permutation Method with the Suwa-Todo Algorithm 

We consider to perform a replica permutation among all M replicas as a generalization 
of a replica exchange in Eq. fll9p : 

X. = {4'"« . ■ ■ • . 4;*'" } X, = {xW"l . . . , } , (22) 

where both i and j are permutation functions and i ^ j- Not only an exchange between two 
replicas, but also a permutation among more than two replicas are allowed in this method. 
Note that the number of all possible combinations between the replicas and temperatures is 
M\. Therefore, the index a of takes a value between 1 and M\. 

In the Metropolis algorithm, the transition probability P {X^ — )■ Xp) for this replica per- 
mutation is also given by Eq. ( 12T]) . The replicas are allowed to transit to non- neighboring 
temperatures, but P (Xq, — > Xp) takes a quite small value for such replica permutation. 
Accordingly, most of such replica permutations are rejected. The number of replica per- 
mutations in which any replica does not transit to non-neighboring temperatures is given 
by 



f] 



^ ^ M-nCn ; (23) 



n=l 



where 

'M' 



^ , for even number of M , 

M (24) 
— 2— , for odd number of M . 

On the other hand, the number of all replica-permutation candidates is Ml — 1 because the 

current combination between the replicas and temperatures is not included. Therefore, the 

probability of trying replica permutations which do not include non-neighboring transitions 

/ [Ml \ 

is I J2n=i M-nCn ) /{Ml — 1). Thus, most of the replica-permutation trials are rejected for 
large M, because 

[f-] M 

J2 M-nCn < J2 ^iCn = 2^' < Ml . (25) 

n=l n=0 
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To avoid this rejection problem, we apply the Suwa-Todo algorithm to the replica per- 
mutations. As in Sec. Ill At we assume that w-r (Xi) is the maximum weight without loss of 
generality. The amount of stochastic flow v {X^ — )■ Xp) is determined by the weight alloca- 
tion in the same way also as in Sec. Ill Al only by replacing the weight Wi to (Xq,). From 
Eqs. ©, (Uni), and (HI]), (X„ ^ Xf^) is given by 

V {Xa = max [0, min [A„^, (X„) + wr (X^) - A^^, wr , (X^)]] , (26) 

where 

Aap = Sa-S(s-i+WKiX^) , (27) 

and 

a 

5, = J]w7r(X^) , So = Smi . (28) 

/3=1 

If wr (X^) (7 7^ 1) is the maximum weight more generally, Eqs. (^7^ and f l25]) are modified 
as follows: 

A^p = Sa- Sfs^i + WK (X^) , (29) 

and 



^ wr (Xp) , for a > 7 



Sa — < 



m7 . (30) 

^ wr (Xp) + '^wk (X/3) , for a < 7 , 

/3=7 /3=1 



5'o = "Sm! • 

A replica-permutation simulation with the Suwa-Todo algorithm is performed as follows: 

Step 1: The label a {a = 1, ■ ■ ■ , M!) of X^ is assigned to all combinations between the replicas 
and temperatures. Table [I] shows an example for M = 3 (three replicas). 

Step 2: For each replica, a canonical MD or MC simulation at the assigned temperature is 
carried out simultaneously and independently for a certain steps. 

Step 3: A replica-permutation trial is performed as follows: First, each weight is obtained 
by Eq. f[T5]) . and the maximum weight wr (X^) is determined. Next, we calculate 
the amount of stochastic flow v (Xq, — )■ X^) in Eq. fl26p and the transition probability 
P (Xq, X/3) = V (Xq, X/3) /wr (Xq,) for /3 = 1, ■ ■ ■ , Ml. Finally, transition from 
State Xq, to State X^ is accepted stochastically with the probability P (X^ — )■ X^). 
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Repeating Step 2 and Step 3, we can carry out the replica permutation MD or MC simulation. 

Figure [2] shows an example of time series of temperatures in RPM. This method realizes 
not only minimization of the rejection ratio but also transitions of the replicas to non- 
neighboring temperatures. 

The number of combinations between the replicas and temperatures increases in propor- 
tion to Ml. For a large number of replicas, we can divide all replicas and temperatures into 
subsets to decrease the number of combinations. Although such a division is not necessary, 
three to eight replicas are appropriate in each subset. As an example, let us consider that 
the total number of the replicas is eight and that they are divided into two subsets which 
have four replicas. In this case, the replica-permutation simulation is performed as follows: 



Step 1: Let us suppose that the temperatures are assigned to the replicas at an initial state as 

^ Replica 1 at Ti ^ 
Replica 2 at T2 

y Replica 8 at Tg y 
They are divided into two subsets: 



^ Replica 1 at Ti ^ 
Replica 2 at T2 
Replica 3 at T3 
y Replica 4 at T4 y 



and 



^ Replica 5 at T5 ^ 
Replica 6 at Tq 
Replica 7 at T7 
y Replica 8 at Tg / 



All combinations between the replicas and temperatures for the former subset and 



those for the latter subset are labeled by and 



[a 



,4:1), respectively. 



Moreover, two more subsets are prepared so that components would differ from those 
of the previous subsets: 



^ Replica 3 at T3 ^ 
Replica 4 at T4 
Replica 5 at T5 
Y Replica 6 at Tq J 



and 



^ Replica 7 at T7 ^ 
Replica 8 at Tg 
Replica 1 at Ti 
y Replica 2 at T2 J 



The labels for the former subset and the latter subset are X^ and X^, respectively. 
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Step 2: For each replica, a canonical MD or MC simulation at the assigned temperature is 
carried out simultaneously and independently for a certain steps. 

Step 3: At an odd number of trail time, replica permutations for and are carried out. 
That is, four replicas are permutated among four corresponding temperatures in each 
subset. At an even number of trial time, replica permutations for X^ and X^ are 
performed. 

Step 4: As a result of the replica permutation for X^, let us suppose that the combination of 
the replicas and temperatures is changed as 

^ Replica 1 at Ti ^ 
Replica 2 at T2 



Replica 3 at T3 
y Replica 4 at T4 y 



^ Replica 1 at Ti ^ 
Replica 3 at T2 
Replica 2 at T3 
Y Replica 4 at T4 y 



(31) 



Due to this permutation, the combination in X^ is automatically changed without a 
replica permutation for this subset: 

^ Replica 2 at T3 ^ 
Replica 4 at T4 
Replica 5 at T5 



^ Replica 3 at T3 ^ 



Replica 4 at T4 
Replica 5 at T5 
y Replica 6 at Tg y 



by Eq. (EI 



(32) 



y Replica 6 at Tg y 



To avoid re-labeling for X^, we perform the next permutation regarding Replica 2 
as Replica 3. Although the combination of the temperatures and replicas may be 
changed by the last replica permutation, we do not need to relabel the combinations 
as in Table [T] by replacing the replicas in this way. 

Repeating Step 2 to Step 4, we continue the replica permutation MD or MC simulation. 



D. Re weighting Techniques 

The results obtained from RPM can be analyzed by the reweighting techniques as in 
REM {sgI, [s^I . Let us suppose that we have carried out a RPM simulation with M replicas 
and M different temperatures (m = 1, ■ ■ ■ , M). 
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For appropriate reaction coordinates and ^2, the canonical probability distribution 
Pt{Ci:^2) at any temperature T can be calculated from 

M 

Pt(6, 6) = E , (33) 



m=l 



and 

e-^" = E^^'"(^i'^2) . (34) 
Ci,6 

Here, = ^ + Sr^, is the integrated autocorrelation time at temperature T^, 
Nm{E] (,1,(^2) is the histogram of the potential energy and the reaction coordinates ,^1 and 
^2 at Tm, and rim is the total number of samples obtained at T^. Note that this probabil- 
ity distribution is not normalized. Equations fl33|) and are solved self-consistently by 
iteration. For biomolecular systems the quantity can safely be set to be a constant in 



the reweighting formulas [37|], and so we set gm = ^ throughout the analyses in the present 
work. These equations can be easily generalized to any reaction coordinates (^1, ^2, ■ ' ' )• 

From the canonical probability distribution Pt(^i,^2) in Eq. (!33|) . the expectation value 
of a physical quantity A at any temperature T is given by 

{A)t = ^ • (35) 

We can also calculate the free energy (or, the potential of mean force) as a function of the 
reaction coordinates ^1 and ^2 at any temperature T from 

FT(ei,6) = -^BTlnPT(6,e2) . (36) 



III. COMPUTATIONAL DETAILS 

A. Asymmetric Double- Well Potential Energy 

In order to verify that RPM gives correct ensembles, we first applied this method to 
a simple system. The system has 100 non-interacting particles with a one- dimensional 
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asymmetric double-well potential energy. The potential energy V (q) at a coordinate q is 
defined by 

V (g) = {{q + if - 1) ((g - 1)^ - 0.9) kcal/(mol ■ A^) . (37) 

To see usefulness of the Suwa-Todo algorithm for the replica-permutation method, we per- 
formed replica-permutation MD simulations with both Metropolis and Suwa-Todo algo- 
rithms. From now on, we will call the replica-permutation method with the Suwa-Todo 
algorithm RPM and that with the Metropolis algorithm M-RPM. The MD versions of the for- 
mer and latter are called RPMD and M-RPMD, respectively. Conventional replica-exchange 
MD (REMD) simulations were also carried out for a comparison purposes. We prepared 40 
different initial conditions (ICs) for each method. The MD simulation from each IC was 
performed for 10.0 ns including equilibration run for 1.0 ns. The time step was taken to be 
1.0 fs. The mass of each particle was 1.0 a.u. Six replicas were used, and temperatures were 
set at 200 K, 235 K, 275 K, 325 K 380 K, and 450 K. The temperatures were controlled by 
the Gaussian constraint method 0, to avoid the problem of non-ergodicity in the Nose- 
Hoover thermostat 6|-l8j for the non-interacting particle system. The trajectory data were 
stored every 10.0 fs. Replica permutations or exchanges were tried every 1.0 ps. 



B. Met-Enkephalin in Vacuum 

To see the usefulness of RPM for biomolecular systems, we next employed a Met- 
enkephalin molecule in vacuum as a test system. The results were compared with those 



obtained from a REMD simulation. The AMBER parm99SB force field 38|, |39| was used. 



The N-terminus and the C-terminus of Met-enkephalin were blocked by the acetyl group 
and the N-methyl group, respectively. Therefore, the amino-acid sequence is Ace-YGGFM- 
Nme. The SHAKE algorithm [4^ was employed to constrain bond lengths with the hydrogen 
atoms during our simulations. The temperature was controlled by the Nose- Hoover thermo- 



stat 



j^. The time step was taken to be 1.0 fs. The initial conformation was an extended 
structure. 

The number of replicas was 12, and temperatures were 200 K, 230 K, 265 K, 300 K, 
340 K, 385 K, 435 K, 490 K, 555 K, 635 K, 720 K, and 820 K. For PRM, we divided the 
12 replicas into two subsets, which had six replicas and six temperatures. The RPMD and 
REMD simulations were performed for 50.0 ns per replica including equilibration run for 1.0 
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ns. The trajectory data were stored every 100 fs. Replica permutations or exchanges were 
tried every 1.0 ps. 



C. C-peptide in Explicit Water 



In order to demonstrate sampling efficiency of RPM for a lager s yste m, we performed 
a RPMD simulation of a C-peptide analog in explicit water solvent 32|]. This peptide is 
known to form an a-helix structure at a lower temperature than 318 K at pH 5.2 32, 



35] . Because the charges at the peptide termini affect helix stability [3l|, we blocked the 
termini of C-peptide by neutral Nme and Ace. Therefore, the amino-acid sequence is Ace- 
AETAAAKFLRAHA-Nme. The histidine residue was protonated to conform our simulation 
conditions to the experimental pH. A REMD simulation was also carried out. The number of 
water molecules was 1800, and two chlorine ions were added as counter ions. The AMBER 
391] was used. The model for the water molecules was the TIP3P rigid- 
emperature was controlled by the Nose-Hoover thermostat 



parm99SB [sM 



body model |41 |. 



SHAKE algorithm 
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. The 



40| was employed to constrain bond lengths with the hydrogen atoms of 



C-peptide and to fix the water molecule structures during our simulations. This system was 
put in a cubic unit cell with the side length of 38.6 A with the periodic boundary conditions. 
The cutoff distance for the Lennard- Jones potential energy was 12.0 A. The electrostatic 
potential energy was calculated by the Ewald method 42(|. The multiple-time-step method 
43I was employed in our MD simulations. For interactions between the C-peptide atoms 
and those between the C-peptide atoms and the water atoms, the time step was taken to 
be 1.0 fs. The time step of 4.0 fs was used for interactions between the water atoms. The 
initial conformation was an extended structure. 

The numbers of replicas in the RPMD and REMD simulations were 24, and temperatures 
were 281 K, 285 K, 289 K, 294 K, 299 K, 304 K, 309 K, 314 K, 320 K, 326 K, 332 K, 338 K, 
344 K, 351 K, 358 K, 365 K, 372 K, 380 K, 388 K, 396 K, 405 K, 414 K, 423 K, and 433 K. 
These replicas were divided into four subsets in RPM. Each subset had six replicas and six 
temperatures. The RPMD and REMD simulations were performed for 40.0 ns per replica 
including equilibration run for 4.0 ns. Namely, the production run of each simulation was 
carried out for 864.0 ns in total. The trajectory data were stored every 400 fs. Trials of 
replica permutations and exchanges were performed every 4.0 ps. 



14 



IV. RESULTS AND DISCUSSION 



A. Asymmetric Double- Well Potential Energy 

We first show that RPM (and M-RPM) gives correct probability distributions for the 
asymmetric double-well potential energy. The probability distributions P (g) in the RPMD, 
M-RPMD, and REMD simulations are presented in Fig. [31 Here, P (g) for each method 
was obtained from the 40 simulations starting from the different ICs. The bin size Aq for 
P (g) was taken to be 0.05 A. The errors were estimated as the standard deviations of the 
40 simulations. To see the accuracy of the simulation results, exact probability distributions 
-Pexact (<?) are also illustrated as the solid lines in the figure. Pcxact (?) at a temperature Tq 
was calculated numerically by 



where Cdw is the normalization constant: Cdw = y dgPc^^ct ■ The integration in 

this equation was computed by the Simpson's method. As shown in Fig. [3l all P (g) show 
good agreement with Pexact (g) at all temperatures. Correct probability distributions are 
obtained by not only REM but also RPM (and M-RPM). 

The transition ratios of the replicas from a temperature to another temperature in each 
method are shown in Fig. ID^a). The transition ratio of the replicas is defined here as a 
probability with which a replica at the temperature is transferred to another temperature. 
RPM has the largest transition ratios at all temperatures among all methods. On the other 
hand, those of M-RPM are extremely small (the values are 0.003-0.004). This is because 
most of the replica-permutation trials were rejected. To realize efficient replica permutations, 
therefore, it is essential to employ the Suwa-Todo algorithm. The total numbers of tunneling 
times of all replicas during the simulations are listed in Table UTI Here, when a replica makes a 
round trip between the lowest and highest temperatures, it is counted as one tunneling. The 
number of the tunneling times is a useful information to see how efficiently the simulation 
samples the temperature space. The value for each method in the table was obtained by 
taking an average of the 40 simulations' results. The number of the tunneling times in RPM 
was 1.7 times larger than that in REM. It means that RPM realizes efficient sampling in 
the temperature space. On the other hand, M-RPM did not have such a large number of 




(38) 
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the tunneling times. This is because its transition ratios were very small. 

To estimate the convergence speed to the exact probability distributions, we examined the 
time series of the average deviations of probability distributions -Ppfrtkie il'j ^) fro^i Pexact il), 
where -Ppfrticie (^'i ^) ^ probability distribution obtained from a single MD simulation of 
particle i from IC j accumulated until time t. The average deviation Ag{t) at the time t is 
defined by 

-j^ ■'Vbin -^particle 

A,(t) = -— -— £^ l^;frticlefe;^)- Pexact (gfc) I , (39) 

iVbiniVlciVparticle ~^ 

where iVbin is the number of the bins, A^ic is the number of ICs, A^particie is the number of 
the particles in a single simulation, and is the coordinate at bin k. 120 bins (A^bin = 120) 
ranging between q^. = —3 and = 3 were taken into account. Figure E] shows Aq(t) 
calculated from each method. Convergence of M-RPM is the slowest due to the low transition 
ratios of the replicas. RPM shows slightly faster convergence than REM at the lowest 
temperature (T = 200 K) although convergence at highest temperature (T = 450 K) is 
almost the same. It is important to increase the sampling efficiency at a low temperature 
because the conformational sampling at a high temperature is originally easier than at a low 
temperature. Therefore, RPM realizes efficient sampling not only in the temperature space 
but also in the coordinate space at a low temperature. It is again because the transition 
ratios of the replicas by RPM are higher than those by REM. 



B. Met-Enkephalin in Vacuum 

The replica-transition ratios and the total numbers of tunneling times for Met-enkephalin 
are shown in Fig. IH^b) and Table [TTl respectively. The transition ratios in RPM are larger 
than those of REM at all temperatures. RPM has also larger tunneling times than REM. 
Thus, efficient sampling in the temperature space was realized by RPM for this biomolecular 
system, too. 

Its local-minimum free-energy structures in vacuum have been reported for various force 
fields although those for the AMBER parm99SB force field have yet to be reported. 
For example, the global-minimum and local-minimum structures for the CHARMM22 force 
field 4J] are shown in Fig. [61 To obtain global-minimum and local-minimum structures in 



the AMBER parm99SB force field and to see the sampling efficiency of RPM, we illustrate 
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free-energy landscapes at 200 K in Fig. [71 The abscissa and ordinate are the root-mean- 
square deviation (RMSD) with respect to the structure in Fig[6]^a) (RMSDi) and that with 
respect to the structure in Fig[H](b) (RMSD2), respectively. Here, the RMSD is defined by 



RMSD = min 



(40) 



where n is the number of the backbone atoms in Met-enkephalin, {q^} are the coordinates 
in the reference conformation, and the minimization is over the rigid translations and rigid 
rotations for the coordinates of the conformation {Qj} with respect to the center of geometry. 
The free-energy landscapes in Figs. El^a) and El^b) were calculated from Eq. fl36|) with the 
reweighting techniques. These landscapes show good agreement with each other, and have 
five local-minimum free-energy states. The five local-minimum states are labeled as A to 
E, as shown in the figures. The free-energy landscapes in Figs. [7](c) and [7](d) were obtained 
from the raw histogram without the reweighting techniques. In the landscape of REM in 
Fig. [Tl^d), States D and E are not observed although these states are observed in RPM. 

In order to discuss sampling efficiency at the lowest temperature more quantitatively, we 
counted the numbers of visiting times in each state, as listed in Table lllli Here, when a 
replica visited a local-minimum state at the lowest temperature after the replica had visited 
another state at the lowest temperature, it is counted as one visit. The errors were estimated 



by the jackknife method 45N47l| in which the production run was divided into 20 segments. 
We regarded the regions presented in Table IIVI as those for the local-minimum states. As 
shown in Table UTTl REM did not visit State D and State E at the lowest temperature. This 
is the reason why these states were not observed in Fig. [7(d). On the other hand, RPM 
visited all states, and the numbers of visiting times in RPM are larger than REM for all 
states. RPM thus samples the conformational space more efficiently than REM at the lowest 
temperature. 

The representative conformations at the local-minimum free-energy states are also shown 
in Fig. [3 The structural features are as follows: (State A) The structure in State A is the 
global-minimum free-energy structure for the AMBER parm99SB force field and similar to 
that for the CHARMM22 force field in Fig. [6]^a). This structure has two hydrogen bonds 
between NH of Gly2 and CO of Phe4 and between CO of Gly2 and NH of Phe4. The 
hydroxy group of the Tyrl side-chain is close to CO of Gly3. (State B) The structure in 
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State B is almost the same as the structure in Fig. Mjo) . Two hydrogen bonds are formed 
between NH of Gly2 and CO of Met5 and between CO of Gly2 and NH of Phe4. However, 
this structure does not have the hydrogen bond between CO of Gly2 and NH of Met5 which 
exists in Fig. MJo). Distance between the hydroxy group of Tyrl and CO of Gly3 is small, 
as in State A. (State C) There are two hydrogen bonds between CO of Tyrl and NH of 
Phe4 and between CO of Tyrl and NH of Met5 in State C. (State D) Two hydrogen bonds 
between CO of Tyrl and NH of Gly3 and between NH of Gly2 and CO of Met5 are formed 
in State D. As for the structures in States C and D, the Tyrl hydroxy group is not close to 
any backbone CO. (State E) The structure in State E has a hydrogen bond between NH of 
Gly2 and CO of Phe4. This structure also has a small distance between the Tyrl hydroxy 
group and CO of Met 5. 



C. C-peptide in Explicit Water 

The transition ratios of the replicas and the total numbers of tunneling times for C- 
peptide in explicit water are shown in Fig. Hl^c) and Table Ull respectively. The transition 
ratios in RPM are larger than those in REM at all temperatures, again. The tunneling times 
of RPM was about 2.1 times larger than that of REM. The time series of the temperatures 
of Replica 1, Replica 9, and Replica 17 are shown in Fig. |H1 Figure E] actually shows more 
frequent tunneling in RPM than in REM. In REM 6 replicas had never made a round trip 
between the lowest and highest temperatures during the simulation. In contrast, all replicas 
had at least one tunneling in RPM. Most of them had more than two tunnelings. Therefore, 
RPM samples the temperature space efficiently than REM. 

It was reported in experiments that C-peptide has a helix structure with salt bridges 



(SBs) between Glu2 and ArglO at a low temperature 33|, |3J|. We obtained helix structures 



which had such SBs in our RPMD and REMD simulations, as in these reports. The lowest 
potential-energy conformation among these helix structures for each simulation is presented 
in Fig. ini Here, we employed the DSSP (define secondary structure of proteins) criteria {3] 
for hydrogen bonds between the side-chains of Glu2 and ArglO and for secondary structures 
of C-peptide. The structure in Fig. |9]^a) obtained from the RPMD simulation has two 
hydrogen bonds between O^ of Glu2 and H^ of ArglO and between O^ of Glu2 and H^ of 
ArglO. The residues from Ala4 to Alall form the a-helix structure. As for the structure 
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from the REMD simulation in Fig. Mjo), the two atoms of Glu2 form the hydrogen 
bonds with the two atoms of ArglO. The a-hehx structure is formed between Ala4 and 
Leu9. Although these structures are slightly different with each other, both RPM and REM 
sampled conformations near the lowest potential-energy helix structure in the other method. 

To see effects of the SBs on the a-helix structures, we calculated probabilities of the 
a- helix structures with the SBs as well as without them. These probabilities at 281 K for 
each residue are shown in Fig. [TOl The probabilities without the SBs in RPM agree well 
with those in REM. The probabilities for residues 4 to 7 in both methods are high regardless 
of the existence of the SBs. This is because their amino-acid sequence is AAAK, and this 
sequence is known for having a- helix structures 49|. In RPM the probabilities with the SBs 
are especially higher than those without the SBs while both probabilities are almost the 
same in REM. We discuss the origin of this difference between RPM and REM later. 

It is considered that the SBs between Glu2 and ArglO stabilize the a-helix structure of 
C-peptide {34 1. In order to investigate the relation between the SBs and the stability of 
the a-helix structure, we calculated a free-energy landscape at 281 K in each method from 
Eq. ( I36l) with the reweighting techniques. These free-energy landscapes are shown in Fig. [HI 
The abscissa is the dihedral-angle distance with respect to a reference a-helix structure. 
Here, a dihedral-angle distance da is defined by 



dn 



-E^K^n, (41) 



nn 



where n is the total number of dihedral angles, Vi is the dihedral angle i, and f° is the 
dihedral angle i of the reference conformation. The distance between two dihedral 

angles is given by 

S{vi,v^) = mmi\vi-v%2n-\vi-v^\) . (42) 

For da, only the backbone-dihedral angles of the residues 4-10 were employed as the elements 
in Eq. ( BTj) . We set the value of to {4>,4') = {—n/S, —tt/3). When the value of da is close 
to 0, therefore, C-peptide has a a- helix structure. The ordinate is distance between of 
Glu2 and of ArglO: Z}(E20e — RIOH^) . Here, this distance is defined to be the smallest 
among the four sets of the distances between two atoms of Glu2 and two H^ atoms of 
ArglO. Six local-minimum free-energy states are observed in both RPM and REM as shown 
in Figs. [TlT a) and (b). We label these local-minimum states as A to F for RPM and A' 
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to F' for REM. The transition states between States A and B and between States A' and 
B', are labelled as G and G', respectively. The a-helix structure with the SBs as in Fig. [9] 
corresponds to a structure in State A or A'. This fact means that the a- helix structure 
with the SBs is a stable structure. On the other hand, the a-helix structure without the 
SBs is not a stable structure because there is no local-minimum states for this structure. 
Therefore, the SBs play an important role in stabilizing the a-helix structure. 

The global- minimum state at 281 K in RPM is State A while that in REM is State 
B'. To see the reason for this difference, we counted the number of visiting times in State 
A for RPM and in State A' for REM at the lowest temperature during the simulations. 
The region from 0.00 to 0.17 for and from 1.0 A to 2.2 A for D{E20, - RIOH^) was 
assigned to State A and State A' here. The number of visiting times for each method is 
listed in Table |Vl Here, we employed two criteria to count the number of visiting times. In 
Criterion 1, when a replica visited in State A (or A') at the lowest temperature after the 
replica had sampled da larger than 0.25 or Z}(E20e — RIOH^) larger than 3.2 A at the lowest 
temperature, it is counted as one visit. In Criterion 2, sampling D{E20^ — RIOH^) larger 
than 3.2 A was not taken into account. Therefore, the number of visiting times increases 
only when a conformation is changed from a non-helical structure to a helix structure with 
the SBs in Criterion 2. In addition to this, breaking and forming of the SBs also increases 
the number of visiting times in Criterion 1. As shown in Table |Vl the numbers of visiting 
times in RPM are much larger than those in REM in both criteria. In Criterion 2, especially, 
the number of visiting times is 1 in REM. Because of such insufficient sampling in State A', 
this state was underestimated in REM. This underestimation caused the small probabilities 
of the a- helix structures with the SBs in Fig. [TOT b). 

The representative conformations for each state in Fig. [TlT a) are shown in Fig. [121 The 
representative conformation for the transition state G is also presented. From these con- 
formations and the free-energy landscape, we clarify a folding pathway from an extended 
structure in State F to the a-helix structure with the SBs in State A. Because the pathways 
in RPM and REM are almost the same, we will discuss that only in RPM: (State F to 
State C via States E and D) The extended structure of C-peptide in State F changes to a 
globular structure as the side-chains of Gly2 and ArglO get close together. Turn structures 
are formed between Gly2 and ArglO by this conformational change, as in States D and E. In 
State C the antiparallel /3-bridge structure is occasionally created between Gly2 and ArglO. 
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(State C to State B) The SBs are formed between Gly2 and ArglO by coming close together. 
Antiparallel /3-bridges between Thr3 and Lys7 or between Leu9 and Hisl2 are occasionally 
observed in State B. (State B to State A via State G) In State G, a short a-helix or 3io-helix 
structure is formed around residues 7 to 10, while maintaining the SBs. By growing this 
helix structure, the longer a-helix structure with the SBs is formed between residues 4 and 
11 as in State A. 

This folding pathway can be divided into two steps. The first step is the "salt-bridge 
formation step", and the second step is the "a-helix formation step". The first step and 
second step correspond to the transition from State F to B and that from B to A, respectively. 
These steps are drawn by the white and red arrows in Figs. [TT] and [T21 We can also see 
in Fig. [TT] that C-peptide rarely takes a folding pathway in which the salt-bridge is formed 
after the a-helix structure formation. 

V. CONCLUSIONS 

We proposed the replica-permutation method (RPM), in which the replicas are allowed 
to transit not only neighboring temperatures but also non-neighboring temperatures. For 
replica-permutation trials in this method, the Suwa-Todo algorithm was employed instead of 
the Metropolis algorithm. This is because most of the permutation trials are rejected in the 
Metropolis algorithm. The Suwa-Todo algorithm had been proposed originally to minimize 
average rejection ratios for state transitions in MCMC. We applied RPM and M-RPM to the 
particles with the double-well potential energy to clarify the usefulness of the Suwa-Todo 
algorithm for the replica permutations. For a comparison purposes, REMD simulations were 
also performed. As a result, RPM realized the most efficient sampling in the temperature 
space while replica permutations were hardly accepted in M-RPM. 

We also applied RPM and REM to Met-enkephalin in vacuum. RPM sampled the tem- 
perature space more efficiently than REM even in the biomolecular system. The five local- 
minimum free-energy states were obtained at 200 K in both methods by using the reweighting 
techniques. In the free-energy landscape estimated from the raw histogram in REM, how- 
ever, two of the five local-minimum states were not observed. This is because REM did not 
sample these two states at the lowest temperature. On the other hand, RPM sampled all 
states even at the lowest temperature. It indicates that RPM realized efficient sampling not 
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only in the temperature space but also in the conformational space. 

Furthermore, the RPMD and REMD simulations were performed for the C-peptide analog 
in explicit water to see the usability of RPM for a larger biomolecular system. RPM showed 
higher sampling efficiency in the temperature space, again. 

It is reported in experiments that C-peptide has the a-helix structure with the SBs 
between Gly2 and ArglO. We observed the a-helix structures in both simulations. We 
also showed that the SBs play an important role in stabilizing the a-helix structure. From 
the free-energy landscape, furthermore, the folding pathway from the extended structure 
to the a-helix structure with the SBs was clarified. This folding pathway consists of the 
two steps. The first step is the "salt-bridge formation step". In this step, the SBs are 
formed by changing its conformation from the extended structure to the globular structures. 
The second step is the "a-helix formation step". The a-helix structure is formed while 
maintaining the SBs. 

We thus revealed that RPM realizes more efficient sampling in the conformational space at 
the low temperature than REM. Furthermore, because the transition ratios of the replicas in 
RPM were larger than those in REM at all temperatures for all systems, larger temperature 
intervals can be taken in RPM. Therefore, the number of replicas can be reduced. 

Although only the results of the MD simulations of RPM were shown in this article, this 
method can be readily applied to the MC method. It is also straightforward to introduce 
this method to the multidimensional REM 50| (also called Hamiltonian REM 5l|) and 
related methods 52|]. We can enhance the sampling efficiency of these methods by replacing 
REM to RPM. 
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TABLE I: Example of an assignment of the labels in RPM with 3 replicas. 



Temperature 


Label a of 


1 


2 


3 


4 


5 


6 


Ti 




Replica 1 


Replica 1 


Replica 2 


Replica 2 


Replica 3 


Replica 3 


T2 




Replica 2 


Replica 3 


Replica 1 


Replica 3 


Replica 1 


Replica 2 


T3 




Replica 3 


Replica 2 


Replica 3 


Replica 1 


Replica 2 


Replica 1 



TABLE IL The total numbers of tunneling times during the simulations. 



Double- Well* Enkephalin C-peptide 
RPM 424 ± 64 459 ±21 58 ± 6 
REM 254 ± 39 357 ± 14 27 ± 4 

M-RPM 11 ±3 

The values were obtained by taking an average of the 40 simulations' results 



TABLE III: The number of visiting times in each state of Met-enkephalin in vacuum. 



Method A B 



C D E 



RPM 97±9 98 ±6 88 ±71±12±1 
REM 86±5 86 ±6 65 ±50±00±0 
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TABLE IV: RMSD ranges for each local-minimum free-energy state. 



State RMSDi (A) RMSD2 (A) 



A 


0.0 - 


0.6 


1.3 - 


1.8 


B 


1.3 - 


1.7 


0.4 - 


1.0 


C 


1.5 - 


1.8 


1.3 - 


1.7 


D 


2.1 - 


2.4 


1.8 - 


2.1 


E 


2.0 - 


2.4 


2.3 - 


2.8 



TABLE V: The numbers of visiting times in State A for RPM and in State A' for REM. 

Method Criterion 1 Criterion 2 
RPM 14 ± 8 5 ± 1 
REM 4 ± 2 1 ± 
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(a) 



14^1 W2 1^3 14^4 

(b) 



W2 
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1/(3-4) 


1/(4-4) 




1/(3-1) 
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FIG. 1: Schematic figures of tlie weight allocation of the (a) Metropolis and (b) Suwa-Todo 
algorithms. Red frame blocks represent rejected flows v{i — >■ i) (i = 1, • • • ,4). 



S 2 

^ ^1 



Only in RPM 



X 












X 







Replica 3 
Replica 1 
Replica 2 
Replica 4 



Time 



FIG. 2: An example of time series of temperatures in RPM. The transitions of replicas in the red 
square frame is not realized in REM. 
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FIG. 3: The probability distributions P (q) of the coordinates q at (a), (c), (e) T = 200 K and 
(b), (d), (f) T = 450 K. These results were obtained from the (a), (b) RPMD simulations, (c), 
(d) REMD simulations, and (e), (f) M-RPMD simulations. The solid lines are the probability 
distributions in Eq. (p8|) . 
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FIG. 4: Transition ratios of the replicas for (a) particles in the double- well potential, (b) Met- 
enkephalin in vacuum, and (c) C-peptide in explicit water. Temperatures are represented as the 
temperature indices. The smallest and the highest indices correspond the lowest and the highest 
temperatures, respectively. 
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FIG. 5: Average deviation Ag{t) in Eq. ^ at (a) T = 200 K and (b) T = 450 K. The solid line, 
the dashed line, and the dotted line show Ag(t) obtained from the RPMD simulations, the REMD 
simulations, and the M-RPMD simulations, respectively. 




FIG. 6: The (a) global-minimum and (b) local-minimum structures of Met-enkephalin in vacuum 
for the CHARMM22 force field. 
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FIG. 7: Free-energy landscapes at T = 200 K obtained by the reweighting techniques from the (a) 
RPMD and (b) REMD simulations and those calculated from the raw histograms obtained by the 
(c) RPMD and (d) REMD simulations. The abscissa is the RMSD with respect to the structure 
in Fig[6l^a). The ordinate is the RMSD with respect to the structure in Fig[6^b). The unit of the 
free-energy landscape is kcal/mol. The labels A to E show the global-minimum and local-minimum 
free-energy states. The representative conformations at these states are also presented. The dotted 
lines denote hydrogen bonds. The figures were drawn by RasMol 53l |. 
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FIG. 8: Time series of the temperatures of Replica 1, Replica 9, and Replica 17. The left- 
hand figures and the right-hand figures are obtained from the RPMD simulation and the REMD 
simulation, respectively. The temperatures are represented as the temperature indices. Index 1 
and 24 correspond the lowest and highest temperatures, respectively. The production runs started 
from the green dashed lines. Red circles indicate the steps at which the replicas reached the highest 
(lowest) temperature after they had visited the lowest (highest) temperature. The red numbers 
present the tunneling times of the replicas. 
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FIG. 9: The lowest potential-energy conformations of C-peptide, which had a-helix structures 
with SBs between Glu2 and ArglO, obtained from the (a) RPMD and (b) REMD simulations. The 
dotted lines denote the hydrogen bonds between the side-chains. The figures were created with 



RasMol 
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FIG. 10: Probabilities of a-helix structures at 281 K for each residue in the (a) RPMD and (b) 
REMD simulations. Solid line and dashed line show the results with and without SBs between 
Glu2 and ArglO, respectively. 
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FIG. 11: Free-energy landscapes at 281 K obtained from the (a) RPMD and (b) REMD simula- 
tions. The abscissa is the dihedral-angle distance with respect to the reference a- helix structure. 
The ordinate is distance between of Glu2 and of ArglO. The unit of the free-energy land- 
scapes is kcal/mol. The labels A to F and A' to F' show the local- minimum free-energy states. 
The labels G and G' are the transition states between States A and B and between States A' and 
B', respectively. 
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FIG. 12: The representative conformations of the local-minimum free-energy states A to F and 
the transition state G obtained by the RPMD simulation. The dotted lines denote the SBs between 
Glu2 and ArglO. Backbone blue color shows turn structures. The figures were drawn by RasMoI 
53|. 
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